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Abstract - The restricted primitive model with nonadditive hard-sphere diameters is shown to 
have interesting and peculiar clustering properties. We report accurate calculations of the cluster 
concentrations. Implementing efficient and ad hoc Monte Carlo algorithms we determine the effect 
of nonadditivity on both the clustering and the gas-liquid binodal. For negative nonadditivity, 
tending to the extreme case of completely overlapping unlike ions, the prevailing clusters are 
made of an even number of particles having zero total charge. For positive nonadditivity, the 
frustrated tendency to segregation of like particles and the reduced space available to the ions 
favors percolating clusters at high densities. 



Ionic soft matter [TJ[2] is a class of conventional con- 
densed soft matter whose interactions are dominated by 
electrostatics crucially affecting its physical properties. 
Among the most popular representatives of such a class of 
materials are natural and synthetic saline environments, 
like aqueous and non-aqueous electrolyte solutions and 
molten salts, including room-temperature ionic liquids, as 
well as a variety of polyelectrolytes and colloidal suspen- 
sions. Equally well known are biological systems of pro- 
teins. 

The simplest fluid modeling an ionic colloidal suspen- 
sion is the Restricted Primitive Model (RPM) [3] a bi- 
nary mixture of uniformly charged Hard-Spheres (HS) for 
which the like-unlike collision diameter between a particle 
of species 1, of diameter cm = a, and a particle of species 
2 of diameter 022 = tr, is equal to the arithmetic mean 
a i2 d = (fii + °22)/2 = a . The two species are of charge 
±q with equal concentrations to ensure charge neutral- 
ity, and the particles move in a medium of fixed dielectric 
constant e. The phase diagram of this model has been 
widely studied both within computer experiments [3HTU] 
and through analytical theories [TTHT5] . 



From these studies emerged how, in the vapor phase 
of this fluid, and thus in the determination of the phase 
diagram, an important role is played by association and 
clustering. In an old paper, |19j one of us studied a more 
general RPM fluid where it is allowed for size nonadditiv- 
ity amongst the particles: the like-unlike collision diameter 
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differs from exf^ by a quantity A 
called the nonadditivity parameter. It was suggested 
through the use of integral equation theories, that such 
a fluid might have a dramatic change of its clustering 
properties. The nonadditivity of the HS diameters docs 
not destroy the simplifying symmetry of the model but 
it introduces modifications of the properties of the pure 
RPM model making it a paradigm for the self-assembly of 
isotropic particles and a challenge to present day theories 
of fluids. There seems to be a lack of literature on this 
subject excepted for Ref. [2"0] . 

In this letter, we reconsider such a model fluid by using 
more direct, highly efficient numerical simulations. In par- 
ticular we analyze the clustering properties outside of the 
gas-liquid coexistence region. As we will see the clustering 
turns out to be greatly affected by the nonadditivity pa- 
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rameter. To the best of our knowledge this is the first time 
that such a model fluid is studied with numerical simula- 
tions. The debate on the importance of clustering in the 
RPM is rejuvenated by studying this new model fluid. 

The model system here considered may be realized 
experimentally through a colloid-star polymer mixture 
where both species are charged |21[|2"2"] and may be rele- 
vant for modeling room temperature ionic liquids |23H26j . 
It is the restricted primitive model (RPM) of nonadditive 
charged hard-spheres (NACHS). The RPM consists of N/2 
uniformly charged hard-spheres of diameter a carrying a 
total charge +q and N/2 uniformly charged hard-spheres 
of the same diameter carrying a total charge —q. The 
spheres are moving in a dielectric continuum of dielectric 
constant e. The interaction between ions of apecies i and 
j a distance r apart is given by 

{+00 r < Uij 
m r>a .. > *,i = 1,2 , (i) 
k B Ter r>a * J 

where (3 = 1 /fc^T with T the absolute temperature and k B 
the Boltzmann's constant, qi the charge of an ion of species 
i. The ions form a mixture of NACHS, i.e. an = 022 = 0" 
and <7i2 = c(l + A), with A > —1 the nonadditivity 
parameter. A thermodynamic state is completely specified 
by the reduced density p* = pa 3 = Na 3 /V, where V is 
the volume containing the fluid, the reduced temperature 
T* = ksTea/q 2 , and the nonadditivity parameter A. 

We used canonical NVT Monte Carlo (MC) simulations 
to study the fluid in a cubic simulation box of volume 
V = L 3 with periodic boundary conditions. The long 
range of the 1/r interaction was accounted for using the 
Ewald method [27] . 

We start from a simple cubic configuration of two crys- 
tals one made of species 1 and one made of species 2 jux- 
taposed. The maximum particle displacement, the same 
along each direction, is determined during the first stage 
of the equilibration run in such a way to ensure an av- 
erage acceptance ratio of 50%. We need around 10 5 MC 
steps (MCS) in order to equilibrate the samples and 10 6 
MCS/particlc for the statistics. 

During the simulation we perform a cluster analysis 
in the vapor phase. After each 100 MCS we determine 
the number N n of clusters made of n particles, so that 
J2 n n ^n = N. We assume [28j[29] that a group of ions 
forms a cluster if the distance r, calculated using periodic 
boundary conditions, between a particle of species i of the 
group and at least one other particle of species j is less 
then some fixed value, i.e. r < Uij + 5a where 5 is a pa- 
rameter 0. In all our simulations we chose S = 0.1 (in 
Ref. [5] a detailed study of the sensitivity of the clustering 
properties on this parameter is carried out for the pure 
RPM fluid). Then we take the average of these numbers 

1 Many different ways of defining a cluster have been proposed 
12.15. 30-321 since the Bjerrum theory 1331 of ionic associations first 
appeared. Our choice corresponds to the geometric one of Gillan 
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(N n ). Here Q n = n(N n )/N gives the probability that a 
particle belongs to a cluster of size n. To establish a cri- 
terion for percolation, we also find clusters without using 
periodic boundary conditions. One of these clusters per- 
colates if, amongst its particles, there are two that do not 
satisfy the cluster condition as a pair, but do satisfy the 
condition if periodic boundary conditions are used. 

In Fig. [T] we simulated the fluid at a temperature 
T* = 0.1 above the critical temperature, T* w 0.05, of the 
pure RPM [6.9. 10 . We sec how, at high density, a positive 
nonadditivity is responsible for a gain of clustering in the 
fluid, which tends to admit percolating clusters also due 
to the fact that a positive nonadditivity pushes the fluid 
at densities closer to the maximum density attainable. It 
is well known that in the neutral nonadditive hard-sphere 
fluid a positive nonadditivity tends to demix the mixture 
at lower densities as A increases |54"l439j . so in our fluid 
we will have a competition between the tendency to demix 
in the neutral nonadditive hard-sphere fluid and the ten- 
dency to cluster in the RPM fluid. At p* = 0.45 both the 
pure RPM and the A = +0.3 have percolating clusters. 
Lowering the density we first reach a state, at p* = 0.3, 
where the negative nonadditivity gives the same cluster- 
ing as the RPM and the positive nonadditivity gives higher 
percolating clustering, then a state, at p* = 0.1, where the 
positive nonadditivity gives the same clustering of RPM 
and the negative nonadditivity a higher one, and finally 
a state, at p* = 0.01, at low densities where a negative 
nonadditivity increases the clustering over the RPM fluid 
and a positive nonadditivity diminishes it. Summarizing, 
for the fixed values of |A| used, we find, in agreement with 
Ref. [19] , that: (a) at high density and positive A we have 
more clustering than in the additive model, (b) at high 
density and negative A we have less clustering than in 
the additive model, (c) at low density and positive A we 
have less clustering than in the additive model, (d) at low 
densities and negative A we have more clustering than in 
the additive model. These points can be explained observ- 
ing that a pair of unlike ions have a higher affinity with 
negative A. Thus, in a bulk phase negative A favors ethe- 
rocoordination. Clusters of a given number of ions tend 
to be smaller when A is negative. As a result, at low den- 
sity (where excluded volume plays a small role) , the extra 
affinity due to negative A enhances cluster formation. By 
contrast, at high densities, the increase in available volume 
from the resulting etherocoordination with negative A has 
an important role, reducing the density-driven imperative 
to form clusters in the negative A case. The same argu- 
ments in reverse explain the behavior of a system with pos- 
itive nonadditivity where now homocoordination at high 
density is favored [T9"] . 

To qualitatively reproduce the curves with non- 
pcrcolating clusters we can use the Tani and Henderson 
clustering analysis [251 HS1 ED] with an inter-cluster con- 
figurational partition function the one of an ideal gas 
of clusters, in reduced units, Z intCT rts (V/a 3 ) Nt , where 
JVf = 53n=i ^ n i s ^ e total number of clusters and we 
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assume to have only clusters made of up to n c particles. 
Then the equations for the equilibrium cluster concentra- 
tions are 



(N n )/N = X n z^/p* , n = l,2,...,n c , (2) 
1 = J^n(N n )/N , (3) 



where 



are the configurational intra-cluster partition 



functions in reduced units with zf 1 *™ = 2 and A(= ap* /2) 
is a Lagrange multiplier to be determined by Eq. ([3]). 
Moreover neglecting the excess internal energy of the clus- 
ters we can approximate z" ltra w (v n /a 3 ) l - n ~ 1 '2 n /n\ where 
v n is the volume of an n— cluster. Assuming further the 
cluster to be in a closed packed configuration we can ap- 
proximate v n sa na 3 /\/2. This simple approximation 
is temperature independent and its usefulness is thereby 
quite limited. 

We checked the size dependence of the curves shown in 
Fig. Q] and saw that when we have no percolating clusters 
the curve was unaffected by a choice of an higher number 
of particles (up to 5000), while the curve changed in pres- 
ence of percolating clusters. In this case we found that a 
common curve is given by (N x )/N with x = n/N 6 [0,1]. 
Then, in order to satisfy the normalization condition, 
1 = £„n((iV„)/A0 « j dxxN 2 {(N x )/N), we must have 
{(N X )/N')/((N X )/N") w {N"/N') 2 for two different sizes 
N' and N". 

In Fig. [2] we show the clustering analysis for the fluid 
with A approaching -1 at T* = 0.1 and p* = 0.45. We 
see how letting A approach —1 this stabilizes small neu- 
trally charged clusters and lowers the degree of dissocia- 
tion a = (Ni)/N. The first stable cluster is the dipole: 
the "overlap" of a positive and a negative sphere. This 
are dipoles of moment qr\2 with 7*12 < c(l + A + 8) which 
may lack a gas-liquid criticality [U]. We clearly have a 
transition from a conducting to an insulating phase as 
A goes from to — 1. We expect that in the limiting 
case of A = — 1 the system we obtain is the neutral HS 
fluid of half the density. This is confirmed by a compar- 
ison of the like radial distribution functions with the one 
of the neutral HS even if the A = — 1 fluid simulation 
rapidly slows down into the frozen configuration of the 
overlapping anions and cations. In order to overcome this 
problem one should alternate single particle moves with 
neutrally charged 2— cluster moves. 

In order to qualitatively reproduce the curve of Fig. [5] 



we need to use Eqs. (HJ)-© with 



v intra ^intra 

where z™f & are the configurational intra-cluster partition 



2 Clearly a proper analysis of the cluster volume would itself re- 
quire a MC simulation 1121 . 



functions of a cluster made of s anions and t cations, 
intra 1 f dr 2 ■ ■ ■ dr s+t 



Sltl 7 0st cr3(.,+t-l) 



(s + t) 



(s+t-1) 



K/K = \ 
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Ab/2 



(X/A' ) min{s ' t} 



(4) 
(5) 



r 2 +\ b /t 



(1+A) 



dr/ 



As/2 



<t(1+A) 



r 2 dr , (6) 



where the configurational integral goes only over the rel- 
ative positions and it covers the region f2 Sjt of s anions 
clusters configuration space, Xb = cr/T* is the Bjcrrum 
length, Roman indeces denote the particle species, Greek 
indeces denote the particle labels, a Roman index with a 
Greek subindex denotes the species of the particle corre- 
sponding to the Greek subindex, and r Ml/ denotes the sepa- 
ration vector between particle p and particle v. Eq. ([5]) is 
justified as follows. Let us call the anions i_ = 1_, . . . , s_ 
and the cations j + = 1 + , . . . , t + . From Eq. (0| follows 

intra 

*'* _ i! 2 



Jn t,t 1=2 
i>j = l i,j= 

-3^zij/ II dr M -tl dr 

J ii+ tin j. i 



°M 1=2 k=l 

t t 

2X B /r i+]+ TT +X B /r i+J _ 



1 



t t 



k=l 



i.j = l 



(7) 



where we approximated e~ Xs / r w 1 which is justified at 
high T* < 1/2(1 + A) or low \ B - Now we observe that 
for example ^i + 2_ = l r i + i_ + i"i_2_| with ri_2_ > & and 



3+As/ri +2 _ 



1. So that for negative nonadditivity we 



can further approximate 



intra 
z t,t 



1 



Jilt . t 1=2 k=l 



i.j = l 



3^)/ n dr M-ii^ + fe- 

Jll *>* 1=2 k=l 



1 

Hp 



i=l 



+A B /r I+i 



oc (2t)< 2 *- 1 ) 
t! 2 



(K/K 



(8) 



where the factor (2t)( 2 * -1 ) takes into account the volume 
of Qt,t- Using the same chain of approximations we reach 
Eq. ((5]). We immediately see how z" 1 1 ra oc K/a 3 becomes 
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Fig. 1: (color online). Clustering properties of the fluid at various values of nonadditivity and density. N n are the number of 
clusters made of n particles. We chose 8 = 0.1. In the MC simulations we used N = 100 particles and a number of MCS= 10 7 . 
The insets allow a = {Ni)/N, the degree of dissociation, to be directly read-off from the graph. 



bigger and bigger as A — > —1 and the same holds for all 
the z™^ & which clearly dominate over all the others z^." tra 
with s t. And this qualitatively explains Fig. [5] 
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Fig. 2: (color online). We show the clustering properties of the 
fluid at T* =0.1 and p* = 0.45 at various values of negative 
nonadditivity approaching — 1. N n are the number of clusters 
made of n particles. We chose S = 0.1. In the MC simulations 
we used ./V = 100 particles and a number of MCS= 5 x 10 7 . 

Sufficiently close to the critical point we determined 
the qualitative change in the behavior of the gas-liquid 
coexistence region by switching on a negative or a pos- 
itive nonadditivity. To this aim we divided the simu- 



lation box into to 3 cubes of side £ = L/m and reg- 
istered the density inside each cell pi = Ni/£ 3 , where 
Mi is the number of particles inside the ith cell so that 

3 

Y^iLi^i — N. Then we calculated the density distribu- 

3 

tion function P m (p) = Y^ILi Pm(pi)/rn 3 [321133], where 
Pm(Pi) is the distribution function for the ith. cell. With 
/ Pm(p) dp = 1. Above the critical temperature the den- 
sity probability distribution function can be described by 
a Gaussian centered at the simulation density whereas be- 
low it becomes bimodal with two peaks one centered at 
the gas density and one at the liquid density. 

We start from an initial configuration of particles of 
random species placed on a simple cubic lattice. We equi- 
librate (melt) the fluid for 10 6 MCS/particle. We then 
sampled the distribution function every 10 MCS. To en- 
hance the efficiency of the determination of the cell density 
distribution, every 10 MCS, we choose the subdivision of 
the simulation box in cells with a random displacement 
r = {r x ,r y ,r z ) with r x ,r yi r z £ [0, L]. And we measured 
the distribution function on runs of 10 6 MCS/particle. 

Choosing m — 2 and N = 100 the results for the fluid at 
a temperature T* = 0.02, p* = 0.2, well within the coex- 
istence region of the pure RPM fluid, and A = 0,±T> 
with V = 10-\10- 2 ,10- 3 are shown in Fig. [3J In 
this case the minimum density that can be registered is 
l/£ 3 = 0.2 x 8/100 = 0.016. We see that the pure RPM 
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fluid shows a density distribution function which has three 
peaks with the first peak, which lies below the minimum 
density, at approximately the low density of the gas phase, 
the second peak at the simulations density p* = 0.2 which 
is due to the fact that the fluid develops surfaces between 
the gas and the liquid phase [H], and the third peak at 
approximately the high density of the liquid phase. We see 
from the figure that increasing T> the middle peak is lost 
first in the positive additive model and then in the nega- 
tive nonadditivc models. Moreover for the biggest V the 
peak of the liquid phase is barely visible. This may be due 
to the fact that one had to choose a proper simulation den- 
sity closer to the density of the liquid [53JI33] ■ ^ e c l earr y 
see how this analysis works like a "microscope" on the de- 
gree of nonadditivity predicting an increase(decrease) of 
the coexistence region for small ncgative(positive) nonad- 
ditivity. This behavior can be explained as follows. Posi- 
tive nonadditivity increases the effective excluded volume 
of ions, thereby reducing the density of the liquid phase, 
and negative nonadditivity does the opposite. 

We believe that our results could be relevant for the in- 
terpretation of experimental work on the phase diagrams 
of room temperature ionic liquids [2S]. In these exper- 
imental systems the liquid-liquid binodals shifted above 
and below the one of the pure RPM are observed depend- 
ing on the kind of solvent. If on the one hand this can 
be ascribed to the different dielectric constant of the sol- 
vent [53] on the other hand it is clear that depending on 
the kind of solvent used the anion-cation contact-pairing 
affinity may vary }45| and thus the different experimental 
ionic liquids considered should be more correctly described 
by comparison not just with the pure RPM but with the 
more realistic primitive model with the addition of either 
a positive or negative size nonadditivity. 

In conclusion, we have performed for the first time a MC 
simulations study of the vapor phase of the RPM with non- 
additive diameters, with particular emphasis on its clus- 
tering properties. A density distribution function analy- 
sis shows how the gas-liquid coexistence region evolves by 
switching on the nonadditivity. A negative nonadditivity 
tends to enlarge the coexistence region while a positive 
one to shrink it. 

From the cluster analysis we where able to distinguish 
between two kind of behaviors for the cluster concentra- 
tions. When we are below the percolation threshold the 
curves for the cluster concentrations as a function of the 
cluster size are independent of the number of particles 
used in the simulation and can be qualitatively explained 
by a simple clustering theory where one approximates the 
clusters to form an ideal gas and the n— cluster as formed 
by n non-interacting particles, for not too small density or 
nonadditivity. When we are above the percolation thresh- 
old the curves depend on the number of particles used in 
the simulation and obey a simple scaling relationship. 

At low density, the negative nonadditivc fluid has higher 
clustering than the pure RPM whereas at high densities 
the positive nonadditivc fluid has a greater degree of clus- 
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Fig. 3: (color online). Cell density distribution function for 
the fluid at T* = 0.02, p* = 0.2 and A = 0, ±V with V = 
10" 3 ,10" 2 ,10 _1 . We used N = 100 and m = 2 with 10 6 
MCS/particle. 



tering. The positive nonadditive fluid is the first one to 
reach the percolating clusters upon an increase of the den- 
sity. This is due to the less space available to the ions, 
for a given density, for positive nonadditivity and to the 
frustrated tendency to segregation of like particles at high 
density. A negative nonadditivity tends to greatly enhance 
the formation of the neutrally charged clusters, starting 
with the dipole, as can be predicted from the simple clus- 
tering theory refined at the intra-cluster level. Traces of 
these features can also be read from an analysis of the 
partial radial distribution function and structure factors, 
which will be presented elsewhere. 

In parallel with the density distribution function analy- 
sis we are currently planning to perform a Gibbs ensemble 
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MC study of the gas-liquid binodal to establish more ac- 
curately the dependence on the nonadditivity parameter. 

We hope that the present study could foster additional 
theoretical and computational studies as well as experi- 
mental realizations of these simple but rich fluids. 

* * * 
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